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ABSTRACT 

The rate of merging of dark-matter halos is an absolutely essential ingredient for 
studies of both structure and galaxy formation. Remarkably, however, our quantita- 
tive understanding of the halo merger rate is still quite limited, and current analytic 
descriptions based upon the extended Press-Schechter formalism are fundamentally 
flawed. We show that a mathematically self-consistent merger rate must be consistent 
with the evolution of the halo abundance in the following sense: The merger rate must, 
when inserted into the Smoluchowski coagulation equation, yield the correct evolution 
of the halo abundance. We then describe a numerical technique to find merger rates 
that are consistent with this evolution. We present results from a preliminary study 
in which we find merger rates that reproduce the evolution of the halo abundance 
according to Press-Schechter for power-law power spectra. We discuss the limitations 
of the current approach and outline the questions that must still be answered before 
we have a fully consistent and correct theory of halo merger rates. 



1 INTRODUCTION 



In current cosmological theory the mass density of the Uni- 
verse is dominated by dark matter. The most successful 
model of structure formation is that based upon the con- 
cept of cold dark matter (CDM). In the CDM hypothe- 
sis dark-matter particles interact only via the gravitational 
force. Since the initial distribution of density perturbations 
in these models has greatest power on small scales, the first 
objects to collapse and form dark-matter halos are of low 
mass. Larger objects form through the merging of these 
smaller sub-units. Consequently, the entire process of galaxy 
formation is thought to proceed in a "bottom- up", hierar- 
chical manner. 

Clearly then, the rate of dark-matter-halo mergers is an 
absolutely crucial ingredient in models of galaxy and large- 
scale-structure formation, from sub-galactic scales to galac- 
tic and galaxy-cluster scales. The Press-Schechter (PS) for- 
malism (Press & Schechter 1974) has long provided a simple, 
intuitive, and surprisingly accurate formula for the distribu- 
tion of halo masses at a given redshift over a large range of 
mass scales and for a vast variety of initial power spectra. 
This formalism states that the number of halos per comov- 
ing volume with masses in the range M — » M + dM is (Press 
& Schechter 1974) 



n(M; t)dM 



2 po S c (t) 
7T M 2 a(M) 
5l{t) 



dlna 



x exp 



2a 2 (M) 



dlnAf 



dM, 



(1) 



where po is the background density and 8 c (t) is the criti- 
cal overdensity for collapse in the spherical-collapse model. 



Here, <r(M) is the root variance of the primordial den- 
sity field in spheres containing mass M on average, ex- 
trapolated to z — using linear theory; it can be deter- 
mined from the primordial power spectrum P(k) that is 
specified by inflation-inspired models for primordial per- 
turbations. For power-law power spectra, P(k) oc k n and 
a(M) oc A'/" (3+n)/6 ; in this case, the Press-Schechter halo 
abundance diverges as n(M) oc M (_9+n)/6 as M — > 0, 
and it is exponentially suppressed above a characteristic 
mass M*(t) determined from the condition S c (t) = er(M*). 
The critical overdensity 5 c (t) is a monotonically decreas- 
ing function of t so that the Press-Schechter distribution 
shifts to larger masses with time, i.e. M» > 0. The fraction 
of cosmological mass in halos of mass M ~ > M + dM is 
Mn(M)dM / po, and, at any given time, most of the mass 
resides in halos with masses M ~ M*. Note that although 
the number of halos diverges as M — » 0, the total mass in 
halos remains finite. 

An elegant paper by Lacey & Cole (1993) — and similar 
work by Bond et al. (1991) and Bower (1991) — extended the 
work of Press and Schechter to determine the rate at which 
halos of a given mass merge with halos of some other mass. 
In addition to providing valuable physical insight, these 
merger rates have extraordinary practical value, having been 
applied to galaxy-formation models, e.g., if galaxy mor- 
phologies are determined by the merger history (Gottlober, 
Klypin & Kravtsov 1991); AGN activity (Wyithe & Loeb 
2003); models for Lyman-break galaxies (Kolatt et al. 1999); 
abundances of binary supermassive black holes (SMBHs) 
(Volonteri, Haardt & Madau 2002); rates for SMBH coales- 
cence (Milosavljevic & Merritt 2001) and the resulting LISA 
event rate (Menou, Haiman & Narayanan 2001; Haehnelt 
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1994); the first stars (Santos, Bromm & Kamionkowski 2002; 
Scannapieco, Schneider & Ferrara 2003); galactic-halo sub- 
structure (Kamionkowski & Liddle 2000; Bullock, Kravtsov 
& Weinberg 2000; Benson et al. 2002; Somerville 2002; Stiff, 
Widrow & Frieman 2001); halo angular momenta (Vivit- 
ska et al. 2002) and concentrations (Wechsler et al. 2002); 
galaxy clustering (Percival et al. 2003) ; particle acceleration 
in clusters (Gabici & Blasi 2003); and formation- redshift 
distributions for galaxies and clusters and thus their distri- 
butions in size, temperature, luminosity, mass, and velocity 
(Verde et al. 2001; Verde, Haiman & Spergel 2002). 

Amazingly enough, however, these merger rates are fun- 
damentally flawed. As we show below, the extended-Prcss- 
Schechter (EPS) formulae for merger rates are mathemati- 
cally self-inconsistent, providing two different results for the 
same merger rate.* The two different merger rates are plot- 
ted in Fig. 1. They are equal for equal-mass mergers but 
increasingly discrepant for larger mass ratios. This ambigu- 
ity will be particularly important for, e.g., understanding 
galactic substructure and for SMBH-merger rates. Even the 
smaller numerical inconsistency for mergers of nearly equal 
mass may be exponentially enhanced during repeated ap- 
plication of the formula while constructing merger trees to 
high redshift. Moreover, the ambiguity calls into question 
the entire formalism, even when the two possibilities seem 
to give similar answers quantitatively^. 

In this paper, we discuss the mathematical requirements 
of a self-consistent theory of halo mergers. As recognised al- 
ready (Silk & White 1978; Sheth 1995; Sheth & Pitman 
1997), the merger process is described by the Smoluchowski 
coagulation equation. This equation simply says that the 
rate at which the abundance of halos of a given mass changes 
is determined by the difference between the rate for creation 
of such halos by mergers of lower-mass halos and the rate 
for destruction of such halos by mergers with other halos. 
The correct expression for the merger rate must be one that 
yields the correct rate of evolution of the halo abundance 
when inserted into the coagulation equation The problem 
is thus to find a merger rate, or "kernel," that is consis- 
tent with the evolution of the halo abundance, either the 
Press-Schechter abundance or one of its more recent N- 
body-inspired variants (Sheth & Tormen 1999; Jenkins et 
al. 2001). 

The apparent simplicity of the mathematical problem, 
which appears in equation (7) below, is in fact quite de- 
ceptive. The Smoluchowski coagulation equation is in fact 
an infinite set of coupled nonlinear differential equations. 
The equation appears in a variety of areas of science — e.g., 
aerosol physics, phase separation in liquid mixtures, poly- 
merization, star-formation theory (Allen & Bastien 1995; 
Silk & Takahashi 1979), planetesimals (Wetherill 1990; 
Malyshkin & Goodman 2001; Lee 2000), chemical engineer- 
ing, biology, and population genetics — so there is a vast but 



* We first discovered this inconsistency in Santos, Bromm & 
Kamionkowski (2002). 

I Extended Press-Schechter theory discusses the correlation of 
peaks in the primordial mass distribution. It is the association 
of such peaks with bound halos, which is not necessarily well- 
defined, that leads to these problems with the derived merger 
rates. 



untidy literature on the subject (although see Leyvraz 2003 
for an illuminating review). It has been studied a little by 
pure and applied mathematicians (Aldous 1999). Still, so- 
lutions to the coagulation equation are poorly understood. 
Furthermore, there is virtually no literature on the problem 
we face: i.e., how to find a merger kernel that, when inserted 
into the coagulation equation, yields the desired halo mass 
distribution and its evolution as a solution. 

In this paper, we present a numerical technique to find 
a merger kernel that yields the correct evolution of a speci- 
fied halo mass distribution. We demonstrate this technique 
by applying it to Press-Schechter distributions for power- 
law power spectra. We regret that at this point we still do 
not have results that can be applied to astrophysical merger 
rates (e.g. those valid for CDM power spectra), although our 
techniques can easily be extended to more realistic cases. 
Moreover, although we have indeed found merger rates that 
are mathematically consistent with the desired halo distri- 
butions, our inversion of the Smoluchowski equation is not 
necessarily unique. As we discuss below, there may be other 
merger kernels that also yield the same halo mass distribu- 
tion. More work must be done to determine how to insure 
that the numerical inversion yields the merger kernel that 
in fact describes the process of mergers from gravitational 
clustering of mass with an initial Gaussian distribution. Nev- 
ertheless, the work presented here may be a first step in this 
direction. 

Below, in §2, we first review the extended Press- 
Schechter calculation and show that it gives mathematically 
inconsistent expressions for the merger rate. In Section 3, 
we then discuss the coagulation equation that needs to be 
solved. Section 4 describes our numerical algorithm for find- 
ing self-consistent merger rates. Section 5 and Figs. 2-8 show 
results of our numerical inversion for a variety of power-law 
power spectra. Section 6 answers some common questions 
about this work, and Section 7 provides some concluding 
remarks and outlines some questions that must still be ad- 
dressed in future work. Appendix A provides an alterna- 
tive formulation of the coagulation equation that makes the 
cancellation of divergences explicit. Appendix B provides, 
for reference, derivations of the two fo the known analytic 
solutions to the coagulation equation. 



2 REVIEW OF THE EXTENDED 

PRESS-SCHECHTER CALCULATION 

The extended Press-Schechter theory (Lacey & Cole 1993; 
Lacey & Cole 1994; Bond et al. 1991; Bower 1991) predicts 
the distribution of masses of progenitor halos for a halo of 
a given mass. By manipulating the equations of this theory 
it is possible to obtain an expression for the merger rate 
of halos of mass Mi with those of mass Mi. Lacey & Cole 
(1993) give the following expression for the probability per 
unit time that a halo of mass M2 will merge with a halo of 
mass Mi = Mi — M2 to make a halo of mass Mf 



d 2 p 
dM 2 dt 



1 



d\na(Mf) 



dlnM 



S c (t) 
a(M f ) 



[l-a 2 (M / )/a 2 (Mi)] 3/2 
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x exp 



a 2 (M f ) a 2 (Mi) 



(2) 



In the extended-Press-Schechter formalism, the abundance 
of halos of mass M is still given by equation (1). The total 
number of mergers between halos of mass M 2 and Mi per 
unit time and per unit volume must therefore be 



R(M 1 ,M 2 ;t) = n(M 2 ;t)- 



d 2 p 



(3) 



dM 2 di' 

Since the merger rate is proportional to both n(M\ ; t) and 
n(M 2 ; t) we define a merger kernel Q(M\, M 2 ;t) (which has 
units of cross section times velocity) through 

R(M u M 2 ;t) = Q(M u M 2 ;t)n(Mv,t)n(M 2 ;t). (4) 

From equations (l)-(4) we can derive 



Q(M u M 2 -t) 



o 2 



pQOi M 2 

1 



din of 



dlnMf 



dhiCT2 



dlnM 2 



(1 - af/aj)^ 



exp 



2 



(5) 



where we have adopted the notation a 2 = a(M 2 ), etc. The 
problem with this merger rate is immediately apparent. 
Clearly R(Mi,M 2 ;t) = R{M 2 ,Mr,t) must always hold; i.e., 
the merger rate must be a symmetric function of its argu- 
ments. However, this is not the case for the above definition 
of Q(Mi,M 2 ;t). This can be seen clearly in Fig. 1 where 
we plot the merger kernels Q(Mi,M 2 ;t) and Q(M 2 ,M 1 ;t) 
for a specific case. Although the differences between the two 
predicted merger rates are small over most of the range of 
M\/M 2 plotted, it is important to note that the discrep- 
ancies may have significant consequences for cosmological 
studies. For example, the merger rate for a mass ratio of 
10~ 5 is uncertain by a factor of 10. This could significantly 
affect the predicted number of dwarf galaxies expected to 
be found within clusters. It is very problematic for predic- 
tions of black-hole mergers detectable with LISA, as many 
of the detectable signals may arise from mergers of very 
unequal masses for which the numerical ambiguity is partic- 
ularly pronounced. Furthermore, even for mergers of nearly 
equal masses, where the two predictions agree more closely, 
repeated application of the merger-rate formula (as occurs 
in the construction of merger trees) will lead to a growing 
divergence between the two predictions. 

The extended-Press-Schechter merger rate is in fact 
only a symmetric function of its arguments in two special 
cases. The first, trivial, case is when Mi = M 2 . The second 
is for the case of a distribution of primordial densities de- 
scribed by a white-noise power spectrum, P(k) oc k n with 
n — 0. In this case the merger kernel reduces to 



Q(M u M 2 ;t) = 



(Mi + M 2 ) 
Po 



(6) 



Our aim is to find kernels Q(Mi, M 2 ;t) which (a) are 
symmetric in their arguments and (b) yield the correct evo- 
lution of the halo distribution h(M) when inserted into the 
coagulation equation. The kernel should also (c) satisfy the 
known statistics of dark-matter-halo merger rates. Below, 
we will illustrate a numerical algorithm that can accomplish 
conditions (a) and (b); we discuss the third condition later. 




10 3 io- 2 

Halo mass ratio 



Figure 1. The dark-matter-halo merger rate computed using the 
extended Press-Schechter formula of Lacey & Cole (1993). Results 
are shown for an n = 1 power-law power spectrum, a halo of 
mass Mi = I0 12 h~ 1 Mq at z = in a Universe with Qq = 0.3, 
Ao = 0.7, <t 8 = 0.9 and T = 0.21. The two lines show the result 
of using the two versions of the merger rate predicted by the 
extended Press-Schechter theory. 



3 BASIC EQUATIONS 

During the process of hierarchical clustering, halos of mass 
M will be created via mergers of pairs of halos of masses 
M' and M — M', and they will be destroyed via mergers 
with halos of any other mass, M'. The rate at which the 
abundance of halos of mass M changes is therefore 



n(M) = 



n(M')n{M - M')Q(M' , M - M')dM' 

-n(M) / n(M')Q(M,M')dM', 
Jo 



(7) 



where the dot denotes a derivative with respect to time, 
the first term on the right-hand side describes halo forma- 
tion, and the second describes halo destruction. Note that 
we have suppressed the explicit dependence of n(M; t) on 
time and the possible dependence of Q(Mi, M 2 ;t) on time 
in equation (7) . This equation is known as the Smoluchowski 
coagulation equation (Smoluchowski 1916). It appears in a 
variety of areas in science in which coagulation processes 
occur. One astrophysical example is the theory of planetes- 
imal growth, in which small objects merge to form larger 
objects. In almost all prior applications, the merger kernel 
Q(M,M';t) is specified by (micro) physical processes (note 
that it has units of cross section times velocity) and the 
equations are then integrated forward from some initial mass 
distribution n(M, t = 0) to determine the mass distribution 
at some later time. 

In our case, however, we know the "answer," the mass 
distribution n(M, t) (either a Press-Schechter distribution, 
an improved version such as the Sheth-Tormen distribution, 
or some other similar distribution determined from simu- 
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lations). We need to find the merger kernel Q(Mi, M2',t) 
that yields the desired evolution of this mass distribu- 
tion when inserted into the Smoluchowski equation. This, 
as far as we know, is an unsolved mathematical prob- 
lem. In principle, some variant of the derivation used by 
Lacey & Cole (1993) that imposes the symmetry constraint 
Q(Mi, M 2 ;t) = Q(M 2 ,Mi;t) might be used to determine 
this merger kernel. It is in fact easy to impose the sym- 
metry constraint with some ansatz, such as an arithmetic 
or geometric mean of the two Lacey-Cole results. However, 
the merger kernel must also be consistent, within the dic- 
tates of the coagulation equation, with the evolution of the 
Press-Schechter halo distribution, and we have not yet been 
able to satisfy this constraint with any analytic approach. 

In the absence of an analytic solution, we attempt to 
find a numerical solution to the problem: i.e., can we numer- 
ically find a merger kernel Q(Mi,M2;t) that when inserted 
into the coagulation equation yields the evolution of the PS 
distribution? The answer, as we discuss below, is yes. To il- 
lustrate the technique, we restrict our attention to the Press- 
Schechter mass distribution because of the simplicity of the 
analytic expressions. Moreover, we restrict our attention to 
power- law power spectra, P(k) oc k n , again for simplicity. 
This has the additional advantage that for the case n = 
we have an analytic solution to the Smoluchowski equation. 
However, our technique can be applied equally well to more 
accurate distributions such as the Sheth-Tormen distribu- 
tion and to other power spectra. 

The Press-Schechter rate of change of halo abundance 
is found by differentiating equation (1) and is given by 



2 po <5 C 



He 2 ' 


8c 




2o\ 


8c 


9 



(8) 



Shifting to a time variable r = — ln<5 c (i) and dimensionless 
mass variable z = M/M*(t) [where M*(t) is the charac- 
teristic nonlinear mass scale defined through the relation 
a(M*) = 8 c (r)] this becomes 



dn 

d7 (Zl) 



2 po -2 8c 

z 1 — exp 



7T M* 



1<T 




\ 5 l-l] 


2 0-f 




9 



(9) 



where a = |dlna/dlnM| = (3 + n)/6. For power-law power 
spectra, a(M) = a-(M*) Z - (3+n)/6 = 5 c z" (3+n)/6 . Therefore, 



dn / 2 po ( _ 9+n)/6 

d7 _ VnW/ 1 



x exp 



1 J3+«)/3 

~2 21 



y (3+n)/3 



The Press-Schechter abundance is simply 



/ \ /2 PO (_9+n)/6 

n(z) = V km: 1 Qexp 



<3+n)/3 



(10) 



(11) 



In these variables, the Smoluchowski equation is 
dn , 



dr 



1 r 

= — / n(z')n(z — z')q{z , z — z')Az 

2 Jo 

- / n{z)n(z')q(z,z')Az , 
Jo 



(12) 



where we have used q to denote the merger kernel in our 
new system of mass and time variables. Our goal is to find 
q(zi, Z2;t). Clearly the solution must be proportional to 



Wn/2M i ,(r)/poa. We therefore choose a system of units 
such that yjir /2M*(r)/ 'pod — 1 to simplify the calculation. 
This choice of units removes the explicit time dependence 
from our merger rate, allowing us to find a function q(zi, 22) 
valid at all times. The time dependence of the merger rate 
is absorbed into a time-dependent system of units instead. 

Before discussing our numerical algorithm, we first 
point out that, for the mass functions we are considering, 
there are divergences in the creation and destruction terms 
in the coagulation equation that cancel. As z — > 0, the mass 
function n(z) oc 2~ 7 , where 7 = (9 — n)/6 is between 2 
and 1 for —3 < n < 3. Thus, if q(z\,Z2) does not vanish as 
one of the arguments approaches zero (which is the case for 
n = 0, and as we argue below, should also be the case more 
generally), then there is a non-integrable singularity at the 
lower and upper limits of the creation term in equation (12), 
and one at the lower limit of the destruction term. These 
divergences cancel, however, if we impose an infinitesimal 
lower mass to the limits of integration. Physically, halos of 
some given (scaled) mass 2 are being created very rapidly by 
merging of halos of infinitesimally smaller mass with halos 
of infinitesimal mass, but they are also being destroyed at 
the same rate by merging with infinitesimal-mass objects. 
Appendix A derives an alternative expression for the coag- 
ulation equation that makes the cancellation explicit. As we 
will see below, these divergences complicate the numerical 
inversion, as the matrix to be inverted will have elements of 
vastly differing magnitudes. 



4 NUMERICAL SOLUTION 

To numerically invert the coagulation equation, we deal with 
a discretized version. We divide z into iV intervals of size 
Az, running from 2 = to z = NAz, labeled by an index i. 
Thus, 2» = iAz, Hi = n(zi), and g i3 = q(zi,Zj). We further 
define j/» = dni/dr. The coagulation equation is then 



1 1-1 



njTii—jqj, 



(13) 



3=1 



Equation (13) can then be re- written in a simple matrix 
form y = K • q. The vector q has iV 2 components, corre- 
sponding to the N x N array q(zi, Zj), while the matrix K 
has dimensions of iV x N 2 . K is the kernel matrix consisting 
of the n's in the above equation. To be explicit, 



y 



jk 



with 



(14) 



(15) 



In practice, we determine K by integrating the terms in 
the Smoluchowski equation over each discrete interval Az, 
with q(zi,Zj) linearly interpolated across this interval. This 
results in cancellation of the divergent terms and is exactly 
correct for the n = case, where q(zi, Zj) is everywhere lin- 
ear. This matrix equation can in principle be solved by a 
suitable inversion method, or by a least-squares minimiza- 
tion to find q. Unfortunately, it is simple to see that the 
equation is ill-determined. We have N linear equations, but 
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N(N + 1) /2 unknowns to determine (since qij is symmetric). 
As such, there will be an infinite number of possible solu- 
tions. We are looking therefore not just for any solution, but 
a sensible one. 



4.1 Regularization conditions 

Since the above equation does not uniquely define q we have 
to apply regularization conditions in order to find a phys- 
ically reasonable q t j. Our goal is to minimize the quantity 
f 2 — jy — K ■ q| 2 . Since this is an ill-determined problem, 
we adopt a regularization condition Ri and instead seek 
to minimize y — K ■ q| 2 and R 2 simultaneously. This reg- 
ularization condition should encapsulate the desired phys- 
ical properties of the solution sought. Specifically, we will 
require that the solution be a smooth function of its argu- 
ments, which seems reasonable for any physically-plausible 
solution, and that 5(21,22) > everywhere, as a negative 
merger rate has no physical meaning. 

Our first regularization condition is therefore 



Ri 



[d 2 q/dz 2 ] 2 + [d 2 q/dz 2 ] 2 dz 1 dz 2 



^,(16) 



where or x is an adjustable parameter. This insures that 
the recovered kernels will be smoothly varying, rather than 
rapidly oscillating. We will return to our second regulariza- 
tion condition shortly. 

The Smoluchowski equation is, of course, linear in 
5(21,22). Since Ri is also linear in 5(21,22) we can use 
straightforward linear algebra to solve for 5(21,22). Since 
the matrix to be inverted can be close to singular in some 
instances, we explore the consequences of optimizing the re- 
sulting solution 5(21,22) using a simple minimization tech- 
nique. Specifically, we aim to minimize the quantity: 



TV- 1 ^ 



+ til, 



(17) 



where Ri from equation (16) is replaced by a suitable dis- 
cretized expression (based upon finite differencing) , and the 
1/(TV — 1) scaling ensures that the relative weight given to 
the two terms in the above is independent of TV. We choose 
Ci = \yi\ to give equal fractional weight to each element in 
the summation. We find that, once the values of TV and Az 
have been chosen to suit the power spectrum under consid- 
eration, the minimization of f 2 leads to almost no further 
improvement in the solution. As such, simple matrix inver- 
sion seems adequate to find 5(21, 22). 

However, in some instances the matrix solution will pro- 
duce solutions for which qij < for certain i and j. These 
are clearly unphysical. As the condition qij > is not lin- 
ear in the qij's we instead apply this condition within our 
minimization routine in such cases. The solution found from 
matrix inversion is used as a starting point for the mini- 
mization. We are then able to find smooth solutions to the 
Smoluchowski equation which are everywhere positive. 

Throughout, we adopt a Rl = 1, which produces smooth 
functions 5(21,22) without limiting our ability to find accu- 
rate solutions to y = K • q. The necessary matrix inversion 
is carried out using LU decomposition. To perform the min- 
imization of f 2 we use a direction-set method (Brent 1973). 
We enforce symmetry of the function 5(21, 22) by re-writing 



the linear equations in terms of the TV(TV + l)/2 independent 
components of q. 



5 PRELIMINARY RESULTS 

Calculations have been performed for power-law power spec- 
tra with n = —2, —1, 0, 1, 2 and 3. For the n = case an 
exact solution is known, 5(21, 22) = (z\ + z 2 ) / V2tt ■ Figures 2 
through 7 show, in their left-hand panels, contour maps of 
the function 5(21,22) recovered by the solution method de- 
scribed above for each value of n. The functions are clearly 
symmetric in their arguments and are all smoothly varying. 
For contrast, grey lines show a geometrically symmetrized 
extended Press-Schechter prediction 5 c ps,sym,G(zi, z 2 ) = 
[5cps(2i, 2 2 )5oPs(^2, «i)] 1/2 , where 5 cP s is the extended 
Press-Schechter merger rate corresponding to equation (5). 
The right-hand panels of these Figures show y{z) predicted 
by Press-Schechter theory, together with the y(z) deter- 
mined from the Smoluchowski equation using q determined 
by the techniques described above and using the arithmeti- 
cally (5oPS,sym,A(2l, 2 2 ) = [5ePS (21 , 2 2 ) + 5oPS (22 , 2l )] /2) 

and geometrically symmetrized extended Press-Schechter 
kernels. (Note that the results for arithmetically and ge- 
ometrically symmetrized extended Press-Schechter kernels 
are indistinguishable in these figures.) In every case we are 
able to find a symmetric, smoothly varying solution which 
solves the Smoluchowski equation. The solutions typically 
differ significantly from the symmetrized extended Press- 
Schechter prediction, which does not solve the Smoluchowski 
equation (except for the specific case of n = 0). 

The results obtained depend upon several factors: 

(i) The relative importance given to solving the Smolu- 
chowski equation and meeting the regularization condition 
when searching for a solution (i.e. the value of a^). 

(ii) The number of elements used in the discretized rep- 
resentation of the Smoluchowski equation and the range of 
2 studied (i.e. TV and Az). (We have adopted the approach 
of making TV as large as possible given constraints on com- 
puting resources, and to make Az as large as possible while 
retaining sufficient resolution to find an accurate solution to 
the Smoluchowski equation.) 

(iii) The nature of the regularization condition(s). 

The first and second factors can be addressed easily given 
sufficient computing resources^ , while the third may require 



t Our current calculations represent the function 5(21,22) by an 
TV X TV grid, with TV = 179. Accounting for the symmetric nature 
of q(zi, 22) this requires TV(TV + l)/2 elements. The matrix to be 
inverted therefore contains TV 2 (TV + l) 2 /4 elements, so memory 
requirements scale as TV 4 . Given that matrix inversion is an TV 3 
process (or, at best, jV 2,807 ) the time required to find a solution 
scales as TV 6 . Thus, merely doubling TV increases memory require- 
ments by a factor of 16 and time requirements by a factor of 64. 
With TV = 179 memory requirements arc of order 2 Gb, while 
finding a solution takes around 3 X 10 seconds on a 3 GHz work- 
station. The choice TV = 179 was therefore dictated by both being 
able to fit the calculation into the memory of a 4 Gb workstation 
and taking a not unreasonable amount of time to run. Further- 
more, the compiler used has a maximum array size which pre- 
vents us from making the matrix any larger, although this could 
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Zj z 



Figure 2. Left-hand panel: The recovered solution for q(zi, 22) for an n = —2 power spectrum is shown as a contour plot (black lines) and 
by the coloured shading. For comparison, we plot the symmetrized extended-Press-Schechter prediction as grey contours. Contours are 
drawn between the values indicated in the figure at equal logarithmic intervals (thinnest line corresponding to lowest value). Right-hand 
panel: The solution, y(z), for an n = —2 power spectrum to the Smoluchowski equation. The dashed black line shows the exact Press- 
Schechter prediction, while the dotted line shows the result found using the Press-Schechter mass function, along with the numerically 
recovered solution for 9(21,22). (Note that these two lines coincide, giving the appearance of a single dot-dashed line.) The solid line 
shows the result obtained by using the Press-Schechter mass function together with the symmetrized extended Press-Schcchtcr prediction 
for 9(21,22) — we plot results for both arithmetically and geometrically symmetrized 9(21,22), but the two are indistinguishable in these 
figures. Note that we plot the absolute value of y(z) and indicate regions where the function becomes negative by thinner lines. 




03468 10 02468 10 

Zj z 



Figure 3. As Fig. 2 but for n = —I. 



consideration of other constraints that any solution which 
is to be considered physically reasonable must meet. With 
the current maximum value of N = 179 we find that our 



be trivially circumvented by using multiple arrays to represent 
our matrices. 



results change at the 10-20% level if we reduce N to 101, 
while changes in A 2 at fixed N similarly lead to changes 
in 5(21,22) at the 10-20% level. Clearly there is a need to 
increase N further, or find ways to perform the inversion of 
the Smoluchowski equation more robustly. We find that our 
results are quite insensitive to or x providing it is not made 
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Figure 5. As Fig. 2 but for n = 1. 



too large or too small (i.e. there is a wide range of gr 1 over 
which the results do not change significantly) . ft seems that 
the simple condition that the function be smoothly varying 
is sufficient to obtain physically plausible solutions. 

Our results can be described moderately well by the 
following fitting formula: 

q(z!,z 2 ) = — 1= (Al[z! + Z 2 ] + A2Z1Z2) 

V 2ir 

X exp (-A 4 [«iZ 2 ] A3 ) , (18) 

where Ai, A2, A3 and A4 are free parameters which we de- 
termine by fitting to the numerically determined q(zi, Z2). 
Table 1 lists the values of these parameters. The final col- 
umn of that table lists the maximum percentage deviation 



between the fitting formula and the numerical results as 
shown in Figures 2 through 7. These fitting formulae are 
valid over the range of 2 shown in the figures, but cannot be 
guaranteed to hold for z's outside of these ranges. In par- 
ticular, for n > the parameter 0,2 is negative, which will 
result in q(zi, 22) < for certain z\, z 2 , which is clearly un- 
physical. We hope to establish a better fitting formula in 
future work. 

It is clear from the right-hand panels of Figs. 2 through 
7 that our solutions for q(zi, Z2) produce the correct rate of 
change of halo abundance. To confirm this we can consider 
the distribution of halo masses at time r, n(z;r), where 
z = M/M,(t). At time t + St the distribution of halo 
masses, using the same definition of z = M/M»(t) and not 
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Figure 7. As Fig. 2 but for n = 3. 



z = M/M*(t+St), is n(z;r+Sr) = /[<5r] 2 n(/[<5r]z; r) where 
f(x) = exp[— 6x/(3 + n)]. Using our recovered 9(21, £2) we 
can evolve n(z; r) to n(z; r + St) for a small step St (such 
as would be used in a merger tree calculation). We show the 
results of this test in Fig. 8, which demonstrates our ability 
to reproduce the evolution of the Press-Schechter mass func- 
tion using our merger rates. (The figure also demonstrates 
how poorly the symmetrized extended Press-Schechter esti- 
mate of the merging rate does.) 

5.1 Some comments on the numerical results 

The numerical and ePS merger kernels appear to be mono- 
tonically increasing functions of mass for power-law indices 



n < 0. However, for power-law indices n > 0, they appear 
to peak at masses near the characteristic mass scale M*. 
As expected, the agreement between the numerical and ePS 
results is exact for n — 0, and they become increasingly dis- 
crepant as n departs from 0. In particular, our numerically- 
recovered merger kernel disagrees with the ePS merger ker- 
nel, and the disagreement increases as n departs from 0. In 
fact, we attempted to find a numerical merger kernel, de- 
manding that the merger kernel be equal to the ePS merger 
kernel for equal-mass mergers. However, our algorithm was 
unable to find a consistent and smooth merger kernel with 
this constraint. This leads us to believe that the ePS merger 
rate is invalid, even for equal-mass mergers where the two 
ePS results agree. 
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Figure 8. The evolution of the Press-Schcchter mass function over a small time step St (the value of which is given in each panel). We 
plot the ratio of the final to initial mass function as a function of mass, z. The dashed line shows the expectation from the Press-Schechter 
theory, while the dotted line shows the result of applying our recovered merger rates to the initial mass function. The solid lines show 
the results of applying the two symmetrized extended Press-Schechter merger rates to the initial mass function (these two lines are 
indistinguishable almost everywhere in this figure). 
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Table 1. Values for the parameters, f A\ to A4, for the fitting 
formula given in equation (18) as a function of power-spectrum 
index, n. The final column lists the maximum percentage differ- 
ence between the fitting formula and the numerical results, when 
using the values given in the table, over the range of z\ and 22 
shown in Figures 2 through 7. 



n 




A 2 


A :i 


A 4 


Max. % deviation 


-2 


0.502 


0.102 


-0.450 


0.118 


15.0 


-1 


0.809 


0.219 


-0.504 


0.027 


8.2 





1.000 


0.000 


0.000 


0.000 


0.0 


1 


1.239 


-0.496 


1.580 


0.009 


14.7 


2 


1.547 


-1.089 


-0.555 


-0.004 


4.9 


3 


1.997 


-2.000 


-0.561 


-0.005 


9.3 



Finally, we point out that Muralidar & Ramkrishna 
(1986) outline a numerical technique for finding merger ker- 
nels, but only under the assumption that the merger kernel 
is homogeneous; i.e., that the merger kernel can be written 
in the form q(zi, 22) = z™f (22/2:1), where m is a power-law 
exponent. We found that we were able to implement this al- 
gorithm and recover the merger kernel for the n — power 
spectrum (for which the merger kernel is indeed homoge- 
neous). We then speculated that for more general power- 
law power spectra, the merger kernel would also be homo- 
geneous. However, we were unable to get the Muralidar- 
Ramkrishna algorithm to return smooth kernels for these 
power spectra. Indeed, our numerically recovered merger 
kernels are not homogeneous for these more general power 
spectra. 



6 FREQUENTLY ASKED QUESTIONS 

Below we attempt to answer several questions which fre- 
quently arise in connection with this work. 

What about three-body encounters? The assumption 
that halos grow through binary mergers, which is implicit in 
the Smoluchowski equation, could be questioned. Since the 
halo abundance diverges at low masses it is possible that 
a halo will experience an infinite number of mergers per 
unit time with halos of vanishingly small mass. However, as 
shown in Appendix A, the creation and destruction rates 
due to these halos cancel exactly. Therefore, we can trun- 
cate the halo mass function at some low mass and obtain 
a good approximation to the true merging history. Put an- 
other way, if the mass distribution is described by discrete 
objects, as in an N-body simulation, rather than a contin- 
uum mass distribution, then there are no three-body merg- 
ers. In such a simulation, there could easily be three or more 
halos that merge within a finite time interval. However, one 
can always find a sufficiently small time step so that any ap- 
parent three-body merger will in fact be a rapid sequence of 
two-body mergers (somewhat like the triple-alpha reaction 
in stellar nucleosynthesis, which is in fact a rapid series of 
two-body reactions). 

What can we learn from N-body simulations? In princi- 
ple, the merger kernel can be determined from N-body sim- 
ulations. In practice, though, this will be challenging. If we 
have N mass bins, then there will be N(N + l)/2 entries in 
the merger kernel, and we must have a sufficiently large sim- 



ulation to provide a statistically significant number of merg- 
ers in each of these bins. Such a simulation will also require a 
very large dynamic range so that a huge number of halos of 
widely varying masses will be resolved. However, there may 
be a simpler avenue. If we carry out a constrained realization 
of a single large halo (see, for example, Springel et al. 2001), 
then the accretion of smaller halos by this single large halo 
can be used to determine the merger kernel Q(Mi, M2) for 
a single large Mi as a function of M2 for M2 <C M\. Such 
a numerical calculation may be helpful in checking the nu- 
merical inversion we have carried out here, and perhaps in 
shedding light on the asymptotic behaviour of Q(M\, M2) 
for M2 <C M\. We plan to pursue such N-body calculations 
in future work — initial investigations using simulations with 
CDM power spectra suggest that our numerical results may 
better match the N-body merger rates than the predictions 
of extended Press-Schechter theory. 

What about fragmentation? In a realistic N-body simu- 
lation, some of the mass in halos might be ejected (e.g., by an 
analog of the "slingshot" effect). If so, then the coagulation 
equation should include fragmentation terms. We believe, 
however, that the fraction of mass ejected will be small, and 
that the inexorable trend in hierarchical structure forma- 
tion is to larger masses; this is also mathematically realized 
in the Press-Schechter and extended Press-Schechter theo- 
ries. We therefore idealize the situation by neglecting any 
such fragmentation and attempt to find a solution to the 
pure coagulation problem first. 

What about "smooth" accretion? Some authors have 
noted that merger trees assembled with the extended-Press- 
Schechter merger rates do not produce Press-Schechter dis- 
tributions, as they should, and have thus attempted to ac- 
count for this discrepancy by the supposition that halos can 
gain mass by accreting diffuse matter from the smooth inter- 
galactic medium, as well as by mergers with other halos. In 
our formalism, what others would call "smooth accretion" 
is described by mergers with very low-mass halos. In other 
words, every mass element is contained in a halo of some size 
(as in standard Press-Schechter theory), and smooth accre- 
tion is taken into account by mergers with halos of extremely 
low mass. 

Is the numerically determined merger rate unique? 
Given our choice of regularization condition — namely that 
the function be as smooth as possible and everywhere 
positive — we do indeed find a unique solution. However, we 
could imagine other regularization criteria, perhaps in ad- 
dition to smoothness and positivity, which would lead to 
different answers. Not all such criteria may allow a solu- 
tion to be found. For example, we find that requiring our 
recovered merger rate to coincide with the extended Press- 
Schechter merger rate for the case of equal mass mergers pre- 
vents us from finding a merger rate which accurately solves 
the Smoluchowski equation. 

Can the merger kernel be determined by looking at the 
properties (e.g., size) of the physical halos? In most applica- 
tions of the coagulation equation (e.g., in planetesimal the- 
ory), the merger kernel is determined by the physical cross 
section and relative velocity of the merging objects; i.e., by 
the "micro" -physics. In our case, however, such a kinetic- 
theory description will not apply, as the system of halos is 
"cold" and halos are more likely to merge with those that 
form nearby, rather than a halo drawn at random from the 
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field. In other words, the merger history is determined by 
the clustering properties of the halos (as encapsulated by 
the power spectrum) rather than by any intrinsic property, 
such as their size. 

What about other algorithms for generating merger 
trees? The fact that the ePS binary merger rates, when used 
to construct merger trees, do not yield Press-Schechter dis- 
tributions is well known. There have been algorithms de- 
veloped that mitigate the shortcomings of ePS by correct- 
ing the ePS merger rates with, e.g., three-body mergers or 
"smooth" accretion (e.g. Somerville & Kolatt 1999; Cole et 
al. 2000). There is no unique algorithm for constructing a 
merger tree (Somerville & Kolatt 1999 explore many of the 
possibilities), but some are more successful than others at 
reproducing the distribution of progenitors of dark-matter 
halos as predicted by the extended Press-Schechter theory. 
It should be noted, however, that the merger rates used in 
these merger trees are not necessarily correct, even if the 
algorithm does ultimately reproduce the correct mass func- 
tion. It should be further noted that all of these algorithms 
rely on the extended Press-Schechter merger rates for major 
mergers which, as we have shown, are very likely incorrect. 
As such, these algorithms do not address the problem we 
consider here. In fact, with our self-consistent merger rates, 
merger trees can be constructed without resorting to three- 
body mergers, smooth accretion, or other kluges. 



7 DISCUSSION 

We have described an approach to find physically reasonable 
estimates of dark-matter-halo merger rates. We describe 
techniques for finding numerical solutions to the Smolu- 
chowski coagulation equation which we use to describe the 
hierarchical formation of dark-matter halos. While the ex- 
tended Press-Schechter theory contains an intrinsic inconsis- 
tency in its predictions for halo merger rates (namely that 
the merger rate of halos of mass Mi with those of mass 
Mi is not the same as that of halos of mass M2 with halos 
of mass Mi) our approach is guaranteed to always produce 
self-consistent merger rates. 

We have presented symmetric, smooth solutions to the 
Smoluchowski equation for the case of dark-matter-halo 
growth through gravitational clustering. These solutions can 
now be checked against the results of numerical simulations. 
The same techniques can also be applied to non-power law 
power spectra (although in general the solution will be some- 
what more complex as the merger rate may depend explicitly 
on time for such power spectra and may contain characteris- 
tic mass scales), and also to halo-abundance evolution rates 
derived from fitting functions designed to reproduce the halo 
mass functions found in N-body simulations. 

A number of questions still remain before we have reli- 
able astrophysical merger rates. First and foremost among 
these is how to insure that the numerically recovered solu- 
tion is in fact the solution that correctly describes mergers 
due to gravitational amplification of an initially Gaussian 
distribution of perturbations. The fact that our numerical 
results are highly insensitive to the details of our smooth- 
ing constraint leads us to believe that our algorithm is in- 
deed converging to a unique physical result, but we cannot 
demonstrate this rigorously. Perhaps future analytic work 



inspired by extended-Press-Schechter-like considerations or 
N-body simulations may provide insight into constraints, 
boundary conditions, and/or asymptotic limits that will al- 
low us to confidently zero in on the correct solution. Beyond 
that, the technique will then need to be applied to obtain a 
merger kernel for a CDM power spectrum, rather than the 
scale-invariant power spectra we have considered here. And 
finally, the merger kernel will ultimately have to be consis- 
tent with a Sheth-Tormen or Jenkins et al. mass function, 
or whatever other mass function is determined to be the 
"correct" one. With these self-consistent kernels, the well- 
known problems with ePS formation-redshift distributions 
(i.e., that they are not positive definite) will be corrected, 
and merger trees will be easily constructed. We plan to ad- 
dress these issues in future work. 

The hierarchical formation of dark-matter halos is a fun- 
damental component of a large fraction of current studies in 
the fields of cosmology and galaxy formation. A mathemat- 
ically self-consistent and quantitatively accurate knowledge 
of halo merger rates is therefore extremely valuable. We be- 
lieve that the techniques developed in this work may provide 
a first step for developing such knowledge. 
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APPENDIX A: ALTERNATIVE 
FORMULATION OF THE COAGULATION 
EQUATION 

Here we provide an alternative form for the coagulation 
equation that demonstrates rigorously the cancellation be- 
tween the divergences in the creation and destruction terms. 
This form of the equation is for the time evolution for the 
fraction of the total mass T{< M) contained in halos of 
mass less than M. This time evolution can be written, 



AF(< M) 
At 



2 Po Jo 



AM' M' 



x / dMm(Mi)n(M'-Mi)Q(Mi,M'-Mi) 

.A) 

i r M 

/ AM'M'n(M') 

Po Jo 

poo 

x / dMm(Mi)Q(M',Mi). (Al) 

.A) 



The first step is to rewrite the first term on the right-hand 
side as 



1 p M p oo p oc 

- / AM' I dMi / AM 2 M'5(Mi + M 2 - M 1 ) 

2 Jo Jo Jo 

xn(M 1 )n(M 2 )Q(M 1 ,M 2 ). (A2) 

We then carry out the M' integral which then changes M' 
in the integrand to Mi + M 2 , and changes the upper limits 



to the Mi and M 2 integrals from oo to M and M — Mi, 
respectively; i.e., this expression is then 

rM pM—M\ 

AM 2 (Mi + M 2 ) 



H du 'l 



xn{M 1 )n(M 2 )Q(M 1 ,M 2 



(A3) 



Since Mi and M 2 enter symmetrically in this expression, it 
can be replaced by 

dMi / dM 2 Min(Mi)n(M 2 )<2(Mi,M 2 ). (A4) 

'o Jo 

We then make the replacements M' — > Mi and Mi — > M 2 
in the dummy variables in the second term in equation (Al) 
and find that the first and second terms are identical, except 
for in the limits for the second integral. Combining the two 
expressions, we find 

f M 

dMiMm(Mi) 



dF(< M) 



dt 



- -J 
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< I AM 2 n(M 2 )Q(M u M 2 ). (A5) 

' M-Mi 

It is then simple to see that if n(M) ~ M~ 7 with 7 < 2 
and Q ~ constant as M — > 0, then the divergence in the 
integrand is integrable. 



APPENDIX B: ANALYTIC SOLUTIONS TO 
THE SMOLUCHOWSKI EQUATION 

In this Appendix we provide an elementary approach to 
two known analytic solutions to the Smoluchowski equation. 
These solutions are not new (see Leyvraz 2003 for a compre- 
hensive review), but we include them here in the hope that 
they may provide some insight into more general solutions. 
We consider the models q^ ^constant, and q^ = i + j. The 
first leads to an exponential mass function for large times 
while the second gives the Press-Schechter mass function for 
an n — power spectrum. 

The Smoluchowski equation is 



hi = - y]qi,iriini, 

i=l 

1 



3-1 00 



qijniUj, 



(Bl) 



where qij corresponds to the discretized coagulation kernel. 

We will consider solutions which satisfy the initial con- 
ditions 



ni = 1, n r — r / 1. 

It will also be useful to define 



(B2) 



M P «=]Tj%«; (B3) 

i.e., the various moments of the number distribution, as well 
as the following generating function, 



(B4) 



3=1 
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We now consider each case individually. 

Constant kernel: q t j = const. Inserting the generating 
function into equation (Bl) gives 

dF F 2 

^ = --FM (t). (B5) 

The initial conditions imply that F(z,0) = cxp(— z), and 
setting z = in equation (B5), gives M (t) = 2/(t + 2). The 
generating function can then be solved for: 

( t + 2)*cxp(Vt(t + 2) - (B6) 
Expanding this gives 

^"^y (B7) 
The asymptotic behaviour of equation (B7) is of the form 

n k ~ t~ 2 g(k/t), t — > co, 
fc = O(t), =4e X p(-277). (B8) 

Sum kernel: qij = i + j. Plugging in for the generating 
function here gives 

dF ldF F dF F 

££. = L?£-Mo(t) - - -MAt). (B9) 
dt 2 dz w 2 dz 2 w v 7 

Setting z = in equation (B9) gives Mo(t) — exp(— 1/2), and 
assuming M\{t) = 1, we have, via the technique of charac- 
teristics (Leyvraz 2003), 

e- t/2 + f)c X p[(e-*/ 2 -l)F]. (BIO) 

The solution is then given using Lagrange's expansion 
(Abramowitz & Stegun 1974), 

n k (t) = ^e-'/ 2 exp[-fc(l - (T^l - e^ 2 )*" 1 . (Bll) 

The large-time behaviour can easily be extracted from 
equation (Bll) as 

7 fc— 1 

n k ~ -^- e - fe e- t/2 , t^oo, fc = 0(l), (B12) 
and, from Stirling's formula, as 
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